Gain Mismatch systematics

Goal

We study the systematics effects induced by gain mismatches for 1 yr Litebird scanning strategy and a fake focal plane with just 2 detector pairs, (total 4 detectors). Also, I 've run a couple of tests simulating South Pole data with Atmospheric sims ( -> Reijo ).

  1. The gain is allowed to drift (I used scramblegains routine ) only in one of the two detectors ( namely $detB$) as a random Gaussian fluctations around 0 and with a 5% width during a fraction of precession period ($5, 10, 96$ min).
  2. To mitigate the systematics due to the gain drifts (i.e. T-> P leakage) we estimate the gain for both $detA$ and $detB$ detector, by regressing with a template ( either Dipole or CMB). The biased coefficient $\hat{ g}$ is related to the true one $g$ by: $$\hat{ g}_i = g_i \frac{1}{1+\sigma^2_{uncorr}/\sigma^2_{corr} }$$ where :
    • $g_i$ is the slope parameter obtained by regressing $d_{t}^i$, TOD of detector $i$, against the signal only template $s_t$,
    • $\sigma_{uncorr} $is the variance of uncorrelated component of the two timestreams, defined as $\sigma^2_{i, uncorr} =Var(d_{t}^i −s_t)/2$,
    • $\sigma_{corr} $, is the variance of correlated one $\sigma^2 _{i,corr} =Var(d_{t}^i) −\sigma^2 _{uncorr} $ .
  1. We then multiply $d_t^B$ by the factor , $\hat{g}_{ A} /\hat{g}_{ B}$

We refer to the TODs with gaindrifts injected as pre-correction, whereas the ones corrected with the methodology described above as post-correction.

Simulation Settings

Simulate 1 yr of observations with LB scanning strategy parameters, α=45, β=50, spinrate=0.1rpm, precperiod=96.5min Focal plane encodes 4 detectors at 140 GHz , FWHM= 36 arcmin, NET= 38.6 μK√s, sample rate= 40 Hz and not spinning HWP.

Choosing the template

we firstly want to quantify how the gain correction procedure is affected by the template chosen to perform the regression . We thus consider signal only simulations encoding as a template for the regression either the Orbital Dipole or the CMB Temperature.

We assess the T->P leakage both at the map and i power spectrum level.

In the case of signal only sims, there is no need to destripe. We therefore looked at binned maps.

In [386]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[386]:
In [405]:
import healpy as hp 
import pylab as pl 




inputclfile ='/Users/peppe/work/satellite_sims/sigonly/inputcl.fits' 
try :
    cl = hp.read_cl(filename=inputclfile)
except IOError:
    cl=hp.anafast(cmbmaps, lmax=lmax)
    hp.write_cl(cl=cl, filename='inputclfile' ) 
hitmap=hp.read_map('/Users/peppe/work/satellite_sims/madam_hmap.fits', verbose=False)
hp.mollview(pl.log10( hitmap), min=3.5, max=4, unit='Log (hitcounts)' ,title='')
cl = pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r0_cl_lensed.dat')
lmax=3*hp.get_nside(cmbmaps[0])
l=pl.arange(2, lmax+1)
Nl =lambda w: (pl.pi /180. /60. * w)**2
wp=7.6 # uK arcmin 
Nl0= Nl(wp)

Regressing w/ Dipole

Let's have a look at the input maps

NB: there is no polarization signal.

In [15]:
dipolemap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/dipole/madam_bmap.fits', 
                       verbose=False)

sigonlymap = hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/dipole_CMB_T/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,7))
hp.mollview(sigonlymap [0], unit='K ', title='Dipole T  ', sub=221 , fig=fig )
cmbT =sigonlymap [0] - dipolemap[0]
hp.mollview(cmbT , unit='K ', title='CMB T  ', sub=222 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q  ', sub=223 , fig=fig  )
hp.mollview(sigonlymap [2], unit='K ', title='U  ', sub=224 , fig=fig  )

Injecting gain drifts ( every 10 min )

In [98]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0], unit='K ', title=r'$  T$  ', sub=131 , fig=fig   )
hp.mollview(driftmap [1], unit='K ', title=r'$  Q$', sub=132 , fig=fig )
hp.mollview(driftmap [2], unit='K ', title=r'$U$', sub=133 , fig=fig  )
cldriftdip=hp.anafast(driftmap , lmax=lmax)

Caption: We get similar non-null signal on Q, and U maps: that's the expected T-> P leakage from gain drifts. In T map the leakage is $\sim 2$ orders of magnitude above the input signal (CMB T + Dipole ). A ring-like pattern is observed in Q and U maps and seems to fall exactly where the dipole changes sign.

Correcting for gain drifts

We now estimate the gains by regressing against the Dipole template.

In [95]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview((driftcorrmap [0]  )    ,unit='K', title=r'$  T  $  ', sub=131 , 
            fig=fig  )
hp.mollview(driftcorrmap [1] , unit='K',title=r'$  Q  $', sub=132 , fig=fig            )
hp.mollview(driftcorrmap [2] ,unit='K', title=r'$  U   $', sub=133 , fig=fig          )
cldriftcorrdip=hp.anafast(driftcorrmap  , lmax=lmax)

Caption: Notice that Q and U maps are a factor ~10 times smaller than the pre-corr. ones (see above ). They present both a dipole-like structures and a couple of scanning strategy features (the two horizontal stripes along the Ecliptic plane and in the Ecliptic poles) . The latter are the regions where pixels are observed $\sim 3 $ times more wrt the rest of the sky, as shown in the hitmap above.

In [371]:
fig=pl.figure(figsize=(16, 5))

pl.subplot(121)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[ 2:lmax+1,2]   ,'k', label='CMB')

pl.loglog(l, cldriftdip[1][2:] *1e12, label='pre-correction')
pl.loglog(l, cldriftcorrdip[1][2:] *1e12, label='post-correction')
pl.legend(loc='upper right')
pl.xlim([2,lmax])

 
pl.subplot(122)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k')
pl.loglog(l, cldriftdip[2][2:]  *1e12)
pl.loglog(l, cldriftcorrdip[2][2:]  *1e12)
pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e0])
foo=pl.xlim([2,lmax])

Caption: Solid blue (orange) are power spectra of pre(post)-gain-correction residual maps. Since no noise has been coadded to the TODs so far, the plots above represent the level of accuracy our estimator might achieve whith the Dipole as a template. As a reference, we plot the CMB spectra from $\Lambda CDM$. We also include in BB plot the solid red line resembling the Litebird whitenoise level of $7.6 \mu $K arcmin, after FG subtraction.

Pros: By comparing the blue with the orange solid lines should highlights the effects of gain correction. As expected they are small for the temperature maps, but in terms of EE and BB power the correction is very effective: it may correct the T->P leakage as much as 3 orders of magnitude in the $\sim$ degree scales.

Cons: the correction doesn't work at $\ell<10$ scales. The reason is related to the presence of Dipole: in facts, it introduces $\ell \geq 2 $ bias in the CMB that needs to be removed by pairing a dipole-fitting code with comp.sep. algorithms (see Litebird CDR, sect. 4.6.1).

Regressing w/o Dipole

We keep the same simulation setup as before, we just changed the input maps encoding now only CMB T (assuming Dipole was correctly subtracted).

In [65]:
sigonlymap = hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/CMB_T/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(12,8))
hp.mollview(sigonlymap[0]  , unit='K ', title=' T  ', sub=131 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q  ', sub=132 , fig=fig  )
hp.mollview(sigonlymap [2], unit='K ', title='U  ', sub=133, fig=fig )

Injecting gain drifts

In [102]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/wodipole/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0]   , unit='K ', title=r'$  T$  ', sub=131 , fig=fig   )
hp.mollview(driftmap [1], unit='K ', title=r'$  Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftmap [2], unit='K ', title=r'$  U$', sub=133 , fig=fig , norm='hist' )
cldriftsT=hp.anafast( driftmap , lmax=lmax)

Caption: The pre-correction maps look as noise fluctuating around $10^{-6} $ K. Note that we've shrunk the colorscale to highlight differences .

We now apply the gain correction and output the maps.

In [103]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/wodipole/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0] , unit='K ', title=r'$  T$  ', sub=131 , fig=fig  )
hp.mollview(driftcorrmap [1] , unit='K ', title=r'$ Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [2], unit='K ', title=r'$ U$', sub=133 , fig=fig, norm='hist' )
cldriftcorrT=hp.anafast(driftcorrmap  , lmax=lmax)

Caption: Post-correction residual maps are one order of magnitude smaller than the pre-correction ones, indicating that the methodology as a clear effect at the map level. Also, we achieve 10 times lower residuals by regressing with CMB T maps than in the Dipole case. A T-like pattern is visible both in Q and U maps .

In [372]:
 
fig=pl.figure(figsize=(16, 5))
pl.subplot(121)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[ 2:lmax+1,2] ,'k', label='CMB')
pl.loglog(l, cldriftsT[1][2:]   *1e12, label='pre-correction')
pl.loglog(l, cldriftcorrT[1][2:] *1e12, label='post-correction')
pl.legend(loc='upper right')
pl.xlim([2,500])
pl.ylim([1e-8,1e-2])

 
pl.subplot(122)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k')
pl.loglog(l, cldriftsT[2][2:]   *1e12)
pl.loglog(l, cldriftcorrT[2][2:]  *1e12)
pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e-2])
foo=pl.xlim([2,500])

Caption: Power spectra pre and post correction in the case of T only sims. Notice that the systematic constamination is lower as compared to the one with Dipole, though even large scale B-modes ($\ell <10 $) of post-correction maps are as much as comparable as the expected level of lensing B-modes at those angular scales. However, we are safely below the LB white noise level after Foreground subtraction (solid red line).

In [438]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[438]:

Regressing with T and E scalar anisotropies (no lensing B-modes)

We include scalar $E$ -modes polarization in the input maps. Also, TE power spectrum is no more null.

In [183]:
sigonlymap = hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/scalar_only/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(14,12))
hp.mollview(sigonlymap[0]  , unit='K ', title=' T  ', sub=131 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q  ', sub=132 , fig=fig   )
hp.mollview(sigonlymap [2], unit='K ', title='U  ', sub=133 , fig=fig )
clS=hp.anafast(sigonlymap    , lmax=lmax)

Caption: TQU maps of scalar anisotropies (no lensing, no tensors).

We scramble gains every 10 minutes

In [113]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/scalar_only/10min/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0]  , unit='K ', title=r'$  T$  ', sub=131 , fig=fig   )
hp.mollview(driftmap [1] , unit='K ', title=r'$  Q$', sub=132 , fig=fig )
hp.mollview(driftmap [2] , unit='K ', title=r'$  U$', sub=133 , fig=fig )
cldriftS=hp.anafast(driftmap    , lmax=lmax)

Caption: pre-correction TQU maps w/ 10 min gain drifts .

In [114]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/scalar_only/10min/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0]  , unit='K ', title=r'$  T$  ', sub=131 , fig=fig  )
hp.mollview(driftcorrmap [1] , unit='K ', title=r'$  Q$', sub=132 , fig=fig )
hp.mollview(driftcorrmap [2] , unit='K ', title=r'$ U$', sub=133 , fig=fig )
cldriftcorrS=hp.anafast(driftcorrmap  , lmax=lmax)

Caption: post-correction TQU maps.

We now plot the EE, BB and TE power spectra.

In [373]:
bl=hp.gauss_beam(lmax=lmax, fwhm=pl.deg2rad(36./60))[2:]
fig=pl.figure(figsize=(16, 5))
pl.subplot(121)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[ 2:lmax+1,2]   ,'k', label='CMB')
pl.loglog(l, cldriftS[1][2:]/bl**2 *1e12, label='pre-correction')
pl.loglog(l, cldriftcorrS[1][2:]/bl**2 *1e12, label='post-correction')
pl.legend(loc='upper right')
pl.xlim([2,600])
pl.ylim([1e-8,1e0])

pl.subplot(122)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k')
pl.loglog(l, cldriftS[2][2:] /bl**2 *1e12)
pl.loglog(l, cldriftcorrS[2][2:] /bl**2  *1e12)

pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e-2])

foo=pl.xlim([2,600])

fig=pl.figure()
pl.title('TE')
ll=l*(l+1)/pl.pi/2.
pl.ylabel(r'$D_{\ell}  [ \mu K^2 ]$', fontsize=16)

pl.plot(l, cl[ 2:lmax+1,4] *ll ,'k')
pl.plot(l, cldriftS[3][2:]*ll  /bl**2 *1e12)

pl.plot(l, cldriftcorrS[3][2:]*ll /bl**2  *1e12)
pl.ylim([-2e2,2e2])
foo=pl.xlim([2,700])

Caption: Power Spectra of pre and post correction maps encoding scalar anisotropies (temp. and polarization). (left panel) E-mode spectra. For consistency, with CMB theoretical spectrum we applied a 35 arcmin gaussian beam, whose effects are evident at higher multiples (right panel). (Bottom panel) TE power spectra (multiplied by the normalization factor , $D_\ell = \ell (\ell +1) C_\ell/2\pi$. Notice that at $\ell > 600$ the pre-correct. power spectrum deviates from the theoretical one whereas the post-correction doesn't.

In [439]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[439]:

MC Signal only sims

To better understand the features observed above we 've run TOAST sims on 20 different CMB realizations (no tensor and lensing B-modes).

We also estimate the transfer function defined as :

$$ \mathcal{F}_\ell = \frac{C_\ell^X}{C_\ell^{input}} $$ with $X$ label indicating either pre or post correction power spectra.

In [367]:
meanS,stdS= pl.load('/Users/peppe/work/satellite_sims/sigonly/scalar_only/10min/sigonly_sims.npy')
meanpre,stdpre= pl.load('/Users/peppe/work/satellite_sims/sigonly/scalar_only/10min/sigonly_sims_uncorr.npy')
meanpost, stdpost= pl.load('/Users/peppe/work/satellite_sims/sigonly/scalar_only/10min/sigonly_sims_corr.npy')
    
bl=hp.gauss_beam(lmax=lmax, fwhm=pl.deg2rad(36./60))[2:]
fig=pl.figure(figsize=(20, 10))
pl.subplot(231)
pl.title('TT')
pl.ylabel(r'$\langle C_{\ell} \rangle_{MC}  [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, meanS[0][2:] *1e12/bl**2, label='no drifts', color='k')
#pl.fill_between(l, (meanS[1][2:] -stdS[1][2:])*1e12/bl**2 , (meanS[1][2:] + stdS[1][2:])*1e12/bl**2 , alpha=.4, color='k' )
pl.loglog(l, meanpre[0][2:] *1e12/bl**2, label='pre-correction')
pl.loglog(l, meanpost [0][2:]*1e12/bl**2 , label='post-correction')
pl.xlim([2,600])
pl.ylim([1e-2,1e3])
pl.subplot(232)
pl.title('EE')
#pl.ylabel(r'$\langle C_{\ell} \rangle_{MC}  [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, meanS[1][2:] *1e12/bl**2, label='no drifts', color='k')
#pl.fill_between(l, (meanS[1][2:] -stdS[1][2:])*1e12/bl**2 , (meanS[1][2:] + stdS[1][2:])*1e12/bl**2 , alpha=.4, color='k' )
pl.loglog(l, meanpre[1][2:] *1e12/bl**2, label='pre-correction')
pl.loglog(l, meanpost [1][2:]*1e12/bl**2 , label='post-correction')
pl.legend(loc='upper right', fontsize=12)
pl.xlim([2,600])
pl.ylim([1e-5,1e-1])
pl.subplot(233)
pl.title('TE')
ll=l*(l+1)/pl.pi/2.
pl.ylabel(r'$D_{\ell}  [ \mu K^2 ]$', fontsize=16)
pl.plot(l, meanS[3][2:]*ll*1e12/bl**2 , label='no drifts', color='k')
#pl.fill_between(l, (meanS[3][2:] -stdS[3][2:] )*ll*1e12/bl**2, (meanS[3][2:] + stdS[3][2:])*ll*1e12/bl**2, alpha=.4, color='k' )
pl.plot(l, meanpre[3][2:]*ll*1e12/bl**2 , label='pre-correction')
pl.plot(l, meanpost [3][2:] *ll*1e12/bl**2, label='post-correction')
pl.ylim([-1.5e2,1.5e2])
foo=pl.xlim([2,700])

pl.subplot(234)
pl.ylabel(r'$\mathcal{F}_{\ell}   $', fontsize=20)
pl.semilogx(l,meanpre[0][2:]/meanS[0][2:] ,label='pre-correction')
pl.semilogx(l,meanpost[0][2:]/meanS[0][2:], label='post-correction')
pl.ylim([0,2])
foo=pl.xlim([2,700])
pl.xlabel(r'$\ell$', fontsize=16)

pl.subplot(235)
pl.semilogx(l,meanpre[1][2:]/meanS[1][2:] ,label='pre-correction')
pl.semilogx(l,meanpost[1][2:]/meanS[1][2:], label='post-correction')
pl.ylim([0,2])
foo=pl.xlim([2,700])

pl.xlabel(r'$\ell$', fontsize=16)
pl.legend(loc='best', fontsize=12)

pl.subplot(236)


pl.plot( l,meanpre[3][2:]/meanS[3][2:])
pl.plot(l, meanpost[3][2:]/meanS[3][2:])
pl.xlabel(r'$\ell$', fontsize=16)

pl.ylim([-2,4])
foo=pl.xlim([2,700])

Caption: (top) Mean TT (left) EE (center) and TE (right) power spectra from 20 MC sims. (bottom) Transfer functions for TT (left) EE (center) and TE (right), as defined above. Pre and post correct. transfer functions similarly fluctuates around $1$, as expected. In particular, at $\ell > 400 $, the gain drift systematics inject extra power to both EE and TE (also TT,at higher $\ell$). The spikes in TE are due to the denominator $C_\ell^{input}$ approaching to very small values.

In [443]:
import scipy as sp
from scipy import interpolate
from scipy import optimize
from scipy import integrate
from iminuit import  Minuit

def get_delta_r(inputcl, lmax=200, wp = 7.6, fsky=1., n_steps=1024):
        lmin=2
        
        ell=pl.arange(lmin, lmax)
        tenscl =interpolate.interp1d (ell, 
                 pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r1_cl_unlensed.dat')[:,3][lmin:lmax])
        lenscl =interpolate.interp1d (ell, 
                pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r0_cl_lensed.dat')[:,3][lmin:lmax ])
        nl_w=interpolate.interp1d (ell, inputcl[lmin:lmax]   )
        clobs   =lambda r,l :r *tenscl(l) + lenscl(l) +  nl_w(l)
        cltheor =lambda r,l:  r *tenscl(l)+    lenscl(l)  +  nl_w(l)
        logP    =lambda r, r_inp,l: -   0.5*(2.*l +1 )*(cltheor(r_inp,l)/clobs(r,l)  +pl.log(clobs(r,l)) -(2. *l -1)/(2.*l+1) *pl.log(cltheor(r_inp,l)))
        #print "Looking for the maximum of L(r)--> "
        func=lambda x: - pl.sum( logP(x, 0, ell)  )
        m=Minuit(func, x=1e-3 , error_x=0.0001, limit_x=(0.,1e-2), errordef=1e-8 , print_level=0)
        m.migrad()
        rmax= m.values['x']
        L       =lambda r: pl.exp ((logP(r,0,ell)   ).sum() - ( logP(rmax , 0, ell)  ).sum() )
        Lmap= map(L, pl.linspace(0,.5e-2, 512) )/L(rmax)
        a=0
        b=1e-2
        integr_steps=n_steps
        Normalization= integrate.quad(L, a=a, b=b , limit=1000, epsrel=1e-8)[0]
        R=0
        dr= (b-a)/integr_steps

        deltar=a
        it_counter=0
        while R<0.6801 and deltar<b :
                deltar+=dr
                R=  integrate.quad(L, a=a , b=deltar, limit=1000, epsrel=1e-8)[0] / Normalization
                it_counter+=1

        deltar/=pl.sqrt(fsky)
        #print "Delta r= %.2f x 10^-3  found after %d iterations at %g per cent C.L.  "%(  deltar*1000,it_counter, R  )
        return Lmap, deltar

Assessing the contamination in B-modes

In [444]:
l2=pl.arange(2, 501)
Lref,deltaref =get_delta_r(  pl.ones_like(l2)*Nl(8.8)  ,lmax=100    )
Luncorr,deltar1= get_delta_r( meanpre[2][2:lmax]*1e12 , lmax=100 , wp=wp   )
Lcorr,deltar2= get_delta_r(meanpost[2]*1e12 ,lmax=100 , wp=wp   )
In [452]:
fig=pl.figure(figsize=(10 , 5) )

pl.subplot(121)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k', label='lensing')
pl.loglog(l, meanpre[2][2:]*1e12   , label='pre-correction')
pl.fill_between(l,  (meanpre[2][2:] -stdpre[2][2:])*1e12  , (meanpre[2][2:] + stdpre[2][2:])*1e12 , alpha=.4 )
pl.loglog(l, meanpost [2][2:]*1e12 , label='post-correction')
pl.fill_between(l,( meanpost[2][2:] -stdpost[2][2:] )*1e12 , (meanpost[2][2:] + stdpost[2][2:])*1e12 , alpha=.4 )

pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3,label='noise %.1f uK arcmin'%wp)
pl.ylim([1e-8,1e-4])
pl.xlabel(r'$\ell$', fontsize=16)
pl.ylabel(r'$\langle C_{\ell} \rangle_{MC}  [ \mu K^2 ]$', fontsize=16)
pl.legend(loc='best')


foo=pl.xlim([2,600])
pl.subplot(122)
 

pl.plot(np.linspace(0,.5e-2, 512), Luncorr, label=r'$\delta r=%.2f \times 10^{-3} $'%(deltar1*1e3) )
pl.plot(np.linspace(0,.5e-2, 512), Lref, label=''+r'$\delta r= %.2f \times 10^{-3} $'%(deltaref*1e3),
        color='red' ,alpha=.5, linewidth=3)
pl.plot(np.linspace(0,.5e-2, 512), Lcorr, label=r'$\delta r=%.2f \times 10^{-3} $'%(deltar2*1e3) )

pl.xlim([0., .5e-2 ])
pl.ylim([0.,1.1 ])
pl.legend(loc='best')
pl.xlabel(r'$r$', fontsize=15)
pl.ylabel(r'$\mathcal{L}(r)$', fontsize=15)
Out[452]:
<matplotlib.text.Text at 0x12e885cd0>

Caption: (left) BB power spectra from 20 MC sims, (solid) average power spectrum, (shaded area) the standard deviation. As a reference we plot as a thick red line the noise level of LB as expected after FG subtraction (assuming $7.6\, \mu K\, arcmin$) and the power spectrum of lensing B-modes. (right) Likelihood of $r$ to assess the level of contamination due to gain drift systematics in terms of $\delta r$ . Notice that the post correction likelihood (solid orange) is narrower than the reference one (solid red).

Summary

  1. The methodology outlined in this posting seems to better work with a CMB temperature template than with the Orbital Dipole. It's not very intuitive given the fact that the Dipole signal can be as much as 10 times larger than the CMB temperature. However, Dipole leaves large scale residuals in polarization maps as large as the scalar E-modes.
  2. Regressing with a temperature map we found that the transfer function $\mathcal{F}_\ell$ of post correction spectra constantly fluctuates around 1 in a verywide range of multipoles $\ell \in [2, 500]$.
  3. Gain drifts simulated with TOAST can potentially lead to a 46 % increase in $\delta r $.
  4. Applying the post-correction represents an overall benefits, both in terms of $\delta r$ and transfer function.
  5. Simulate w/ a larger number of detectors doesn't change the features presented in this posting (e.g. tried a couple of runs with 8 detectors to have a better attack angle coverage).

Future Outlooks

  • Inject gain drifts with a typical time scale, estimate and correct gain drifts on a different one
  • Simulate w/ a larger number of detectors surely help to test the performance of this methodology on realistic datasets with white and 1/f noise, (signal+ noise power spectra from sims 2 detector pairs are mostly dominated by noise).